rm(list = ls())

## Packages

library(tidyverse)
library(stargazer)
library(lfe)

## Def helper function

str_filter <- function (string, pattern, neg = F) 
{
  if (neg == F) {
    string %>% .[str_detect(., pattern)]
  }
  else {
    string %>% .[!str_detect(., pattern)]
  }
}

## Get data

df <- read_rds("data/Data_Clean.rds") %>% 
  mutate(year = as.numeric(year)) %>% 
  mutate(unit = coalesce(ags, vb_key)) %>% 
  mutate(treat_weakly = ifelse(treat_categ == 'schwach', 1, 0),
         treat_heavy = ifelse(treat_categ %in% c('schwer', 'sehr schwer'), 1, 0)) %>% 
  mutate(treat_weakly_post = treat_weakly * post,
         treat_heavy_post = treat_heavy * post)

## Get salience, do median split 

sal <- read_rds("data/civey_salience_county_aggregated.rds") %>% 
  dplyr::rename(county_id = ags) %>% 
  mutate(across(matches('salience'), 
                list(~ifelse(.<median(., na.rm = T), 1, 0)), 
                .names = "{.col}"))

## Merge to main data

df <- df %>% 
  left_join(sal)

#### Function to generate differenced DFs ####

do_diff <- function(df, y, diff_years = c('2017', '2021')) {
  
  df_out <- df %>% 
    arrange(unit, year) %>% 
    dplyr::select(unit, year, matches('treat'), !!y, county_id,
                  land, matches('salience'), one_of(covars)) %>% 
    filter(year %in% diff_years)
  
  df_out[, 'y'] <- df_out[, y]
  
  ## Is every unit in there twice?
  
  df_out$unit %>% table() %>% unique()
  
  df_out <- df_out %>% 
    group_by(unit) %>% 
    summarise(y_diff = y[2] - y[1],
              treat_weakly = first(treat_weakly),
              treat_heavy = first(treat_heavy),
              county_id = first(county_id),
              land = first(land),
              county_age = county_age[2] - county_age[1], 
              county_unem = county_unem[2] - county_unem[1], 
              county_pop_tot = county_pop_tot[2] - county_pop_tot[1], 
              county_foreign_share = county_foreign_share[2] - county_foreign_share[1], 
              county_abitur_share = county_abitur_share[2] - county_abitur_share[1],
              treat_wounded_or_dead = first(treat_wounded_or_dead)) %>%
    ungroup() %>% 
    dplyr::select(unit, y_diff, matches('treat'), 
                  county_id, land, matches('county')) %>% 
    dplyr::rename(!!y := y_diff) %>% 
    filter(!is.na(!!y)) %>% 
    distinct() %>% 
    left_join(sal)
  
  ## Return 
  
  df_out
  
}  

#### Table 2 - HTE by salience #### 

## Generate differenced DFs

covars <- df %>% 
  colnames() %>% str_filter("county")

diff_df <- do_diff(df, 'greens') %>% 
  mutate(greens = greens * 100) %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) 
  
## Regressions

m1 <- felm(greens ~ treat_weakly*salience_county + treat_heavy*salience_county | 0| 0| county_id,
           data = diff_df)
m2 <- felm(greens ~ treat_either*salience_county | 0| 0| county_id,
           data = diff_df)
m3 <- felm(greens ~ treat_wounded_or_dead*salience_county | 0| 0| county_id,
           data = diff_df)
m4 <- felm(greens ~ treat_weakly*salience_county + treat_heavy*salience_county | land| 0| county_id,
           data = diff_df)
m5 <- felm(greens ~ treat_either*salience_county | land| 0| county_id,
           data = diff_df)
m6 <- felm(greens ~ treat_wounded_or_dead*salience_county | land| 0| county_id,
           data = diff_df)

## ##

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)


stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     #covariate.labels = c('Weakly affected',
                     #                     'Highly affected',
                     #                     'Any wounded or dead'),
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: any',
                                          'Any wounded or dead',
                                          'Flooding: heavily affected',
                                          'Flooding: weakly affected * low salience',
                                          'Flooding: highly affected * low salience',
                                          'Flooding: any * low salience',
                                          'Any wounded or dead * low salience'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)),
                     type = 'text')

#### Table A.6 - HTE by salience with covars #### 

m1 <- felm(greens ~ treat_weakly*salience_county + 
             treat_heavy*salience_county + 
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share| 0| 0| county_id,
           data = diff_df)
m2 <- felm(greens ~ treat_either*salience_county + 
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share| 0| 0| county_id,
           data = diff_df)
m3 <- felm(greens ~ treat_wounded_or_dead*salience_county+ 
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share| 0| 0| county_id,
           data = diff_df)
m4 <- felm(greens ~ treat_weakly*salience_county + treat_heavy*salience_county + 
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share| land| 0| county_id,
           data = diff_df)
m5 <- felm(greens ~ treat_either*salience_county + 
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share| land| 0| county_id,
           data = diff_df)
m6 <- felm(greens ~ treat_wounded_or_dead*salience_county+ 
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share| land| 0| county_id,
           data = diff_df)

## Models

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)

## Table 

stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     #covariate.labels = c('Weakly affected',
                     #                     'Highly affected',
                     #                     'Any wounded or dead'),
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: any',
                                          'Any wounded or dead',
                                          'Flooding: heavily affected',
                                          'Flooding: weakly affected * low salience',
                                          'Flooding: highly affected * low salience',
                                          'Flooding: any * low salience',
                                          'Any wounded or dead * low salience'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)),
                     type = 'text')

